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Abstract 

In this work we present a theoretical study on the propagation of light in heterogeneous sys¬ 
tems with fluctuating optical properties. To understand the consequences of the fluctuations 
we perform numerical calculations with uniform and non uniforms systems using Monte Carlo 
simulations. We consider two distributions to represent a non-uniform medium: delta func¬ 
tion and an exponential negative distributions.The results show that even with finite moments 
distributions, may require a large number of interactions for a convergence towards Gaussian 
statistics. This can be important when estimating the optical properties of thin films. 
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FIG. 1. Schematic view of a non homogeneous 
particle distribution. 
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I. INTRODUCTION 

Inhomogeneous medium is characterized by having variations in refractive index or have 
inclusions of particles with another index. Normally, it is assumed that the media has 
properties which do not vary with average position. However, many of these systems, 
for example, composed of a suspension of particles, tend to form lumps or regions where 
the particle density is higher (by sedimentation, for example). This breaks down with 
the assumption of uniformity of the system (see Figure 1) and involves certain difficulties 
to study. A major consequence is that, depending on the type of fluctuation, random 
process may not be stationary and hence not ergodic. This means that sample averages 
aren’t equivalent to ensemble averages (ensemble). This class of disordered systems are 
known as superdifusive media mn and have raised recent interest in different areas of 
science, in particular in optics |1]. 

Figure illustrates the type of paths we can expect in a medium with particles of many 
sizes. We see that, in regions with small particles, the light follows a zigzag path with 
step short while to encounter a big particles, abruptaly, the paths are so long. This means 
that in the media there is a non-negligible probability that the light suddenly run, a much 
longer than the other. 

The paper is organized in the following way. In the first section we describe in broad 
terms the problem we will address. Section II discusses the stochastic properties of a 
medium where the statistics are Gaussian and non-Gaussian. Section III presents a study 
of the resultant of flight on uniform and disordered media. Finally, the section provides 
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FIG. 2. Illustration of possible optical paths in a non-homogeneous media with fractal charac¬ 
teristics due to the dispersity in the particle sizes [Ij. 

a summary and conclusions of this work. 


II. STATISTICAL FLUCTUATIONS OF THE OPTICAL PROPERTIES 

In simplest theories, which are approximations to the radiative transport equation, inho¬ 
mogeneous medium can be characterized by their interaction coefficients and /r^, and 
also by the anisotropy parameter g. These coefficients represent, respectively, the scat¬ 
tering probability per unit length (/Xs), the probability of absorption per unit length (pa) 
and average scattering angle after an interaction {g =< cosOg >). The total interaction 
coefficient /Xi = /is + l^a and the mean free path is / = l//is. 

In these simple theories it’s has, moreover, that fig and g always appear as /x'^ = /is(l — g), 
which is known as the reduced scattering coefficient. Henceforth, we assume that the 
medium is not absorbing, so that we can /Xq = 0 and /x* = /Xs and also it’s assume that 
g = 0 (isotropic scattering). 

Generally it’s assumed that the parameter /Xg is constant. However, this is not always 
true in experimental systems. Under the assumption of independent scattering of the 
particles, the interaction coefficient can be written as 

fit = pCt, [cm“^] (1) 

where Cg denotes the scattering cross section of particles and p indicates the bulk density, 
he., the number of particles per unit volume. We can see that if there are changes in the 
density, or if the particles have different properties, pg could be a function of position. 
We can see from equation ([^ that pt depends on the optical size of the particles Ct and 
of the density p. We start the study considering fluctuations with respect of one the 
parameter these parameters. Then we will consider the effect of the fluctuations of both 
of them in the total interaction coefficient. 
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A. Uniform systems 


We consider first the case of a non-homogeneous system on a microscopic level bnt that, 
beyond a certain level, not changed in their optical properties. We model the media using 
identical particles with a uniform particle density. The probability of interaction per unit 
length along a line is then constant. 

Let /is the probability of interaction per unit length of a photon in the media. This 
probability can be written as the sum of the probabilities of scattering and absorption: 
/it = /is + /ia- For simplicity, we consider that there is no absorption in the media so that 
/if = /is. We also dehne the mean free path between interactions as I = (/is)~^- 
Let F{s) the probability that a photon, which begins at s = 0, does not scattered on 
the length s. The probability of scattering in a differential length ds is figds, so that the 
probability of no scattering at ds is (1 — figds). We then have 

F{s + ds) = F{s){l — figds), (2) 

so we can write the differential equation 

dF{s) = —F{s)figds. (3) 

The solution, gives the probability that a photon is not scattered in the length s, 

F{s) = exp{-figs), (4) 

so that the probability of scattering at this length is given by 

P{s) = 1 - exp{-figs), (5) 

The probability density function (PDF), pg{s) = dP{s)/ds, which governs the interaction 
of photons with the media may then be written in the form 

Pg{s) = figexp{-figs) = y exp(-s//), (6) 

The moments of the distribution are given by [7] 

(7) 

so that 

(s) = /, < >= 2f, , (8) 

and standard deviation 

= V(s^) - {sy = 1. (9) 
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B. Systems with two types of particles 


In a non-uniform system it’s present variations of the interaction parameter. We call this 
parameter fluctuating u, which its average is /i^. 

The probability density function of displacement in a given region will depend on the 
specihc value of the random variable u takes. From the equation (§ , then the conditional 
PDF can be writen as 

Psisl^) = z^exp(—z/s). (10) 

Denoting by the PDF to scattering coefficients, we can write an expression for the 

new PDF to displacements, 

Ps{s) = / Ps{s\v)p^{v)dv. (11) 


In the case where the scattering coefficient, take just two random values, we can write 
the PDF as, 

p^iy) = a5{y - pi)+ b5{y - pi 2 )', (12) 

where a is the probability of the coefficient pi and b represents the probability of occurring 
P 2 - It’s necessaryy that a -|- 6 = 1 and ps = o.pi -|- bp 2 . 

At this point, it is necessary to mention that the modeled system is not a homogeneous 
mixing of two components (in which case, ps would be equal to the sum pi -|- P 2 ), but a 
system with regions with properties pi and other with properties p 2 - 
It’s convenient to dehne also a parameter a, by the relationship 


a 

/^1 


( 13 ) 


implies that. 


P2 = 


1 — a 
1 — a 


ps 


(14) 


This means that the PDF (12) can be specihed by the parameters a, pi and p 2 or, 
alternatively, by /is, a and a. 

The PDF for the displacements, or flights, can be determined by equations (119, (0 and 


(12), and is given by 


Ps(s) = j ve [a5{v - pi)+ b5{v - p 2 )]dv, 

= apie~^^^ + bp2e~^^^ (15) 

We can verify that this FDP is normalized properly and that the hrst moments are 


, , a b 

{s) = — + 

Pi P 2 


. 2 , 2a 26 


( 16 ) 
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FIG. 3. Comparing the probability density function defined by equations (15) and (19). The 


solid blue curve represents the distribution for flights governed by a variation in v double-delta 
type. The dotted red curve for the negative exponential distribution. 

C. Systems with negative exponential fluctuations 

For a medium with a distribution of scattering coefficients negative exponential type, we 
have that the PDF for v can be written as, 


Vv[v) =/?exp(-^!/), 


( 17 ) 


where /3 = l//rs and ps is the average scattering coefficient of the system. 

As in the previous case, the equations 0, 0 and we have that the PDF for flight 
is given by, 

poo 

Ps{s) = / [uexp{-us)][/3exp{-(]u)]du, 


= l3 u exp[—(s +/3)iy] du. 

Jo 


Evaluating the integral, we have 


Ps{s) = 


(9 


ps 


(s + ,3)'i (l+p,s)2 


( 18 ) 


( 19 ) 


We see that for large s arguments, the probability density behaves as a Lorentzian Levy 


flight [equation (??)]. It is worth mentioning that the PDF dehned by equation (19) has 
no dehnite time, which is characteristic of Levy flights. 

Figure]^ shows the behavior of the density (15) and (19). The curves were scaled indepen¬ 
dently to illustrate the differences. We can see that, although the curves appear similar, 
they have important differences. In particular, it should be noted that the decay of the 
curve corresponding to negative exponential fluctuations is very slow. 
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FIG. 4. Random walk trajectories, (a) Random walk in a uniform media, (b) Levy random 
walk. 

To better understand the consequences of adopting these PDF, in the next section we 
present calculations of random walks using Monte Carlo simulations. 


D. Fluctuations in density and optical size 

A realistic system must consider fluctuations in the optical size Cx and density p. 

III. RANDOM WALKS 

In the context of this study, it is interesting to see the result of random walks with different 
probability density laws described in the previous section, focusing on situations in which 
the number of interactions is large. The types of situations that may occur in a case with 
uniform distribution and one in which they occur Levy flight is illustrated in Figure 5. 
We consider the distance after N number of displacements. For simplicity, we illustrate 
the method by considering a two-dimensional space and we’ll write the total displacement 
in polar coordinates (a, 6*) (see Figure [^: Evaluating the integral, we have 

( 20 ) 

We assume that: 

1. Amplitudes Sk/VN and phases 0^ are statistically independent. 

2. The variables Sk follow the distribution ([^ with moments given by equation ([^. 

3. The phases cpk are uniformly distributed in the interval (—7r,7r). This means that 
the scattering of the particles is isotropic. 


a = = 


Vn 






4^k 
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FIG. 5. Random walk. 


We then have that the x and y components are given by: 


acos9 = SfcCos(d>fc), 

(21a) 

1 ^ 


a sin 0 = ^ Sk sin(0fc), 

(21b) 


and, with onr assnmptions, we hnd that 

(flx) 0 ) 
{ay) = 0 


G) = 2 . 

G) - 


(22a) 

(22b) 


When the nnmber of steps, iV, is very large, the displacement of the photon statistics 
are Ganssian. That is, both Qx and ay follow as Gaussian distributions.In this case, the 
variables have zero mean, second moment 1/2 and they are not correlated. It is then 
circulated Gaussian aletorio process [7]. The joint probability density is then. Evaluating 
the integral, we have 


^y) ^ exp 



(23) 


The length statistics are founded with a probability transformation to express (23) in 
terms of (a, 9), and integrating over phase [7], It is found that 


Pa(a) = — exp 



(24) 


for a > 0. Then the steps follow a a Rayleigh distribution. 

We note, however, that our assumptions are valid, moments of the distribution that 
governs the movement must be hnite. 
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In the next section, we will use Monte Carlo simulations to study three cases, correspond¬ 
ing to the displacement PDF given by equations (§. ( [I5| and ( p!^ , starting with the case 
of the uniform system. 


IV. MONTE CARLO SIMULATIONS 

To explore the convergence to Gaussian statistics for the three PDF considered, we present 
calculations based on the MCML (Monte Carlo Multi Layered) simulation [9] using the 
values for the average properties of the medium shown in Table I. Table dehnes the 
parameters of the hypothetical medium that we study. Assume that the medium is highly 
scatterer (/i^ = 1867 cm“^) and with low absorption (/Xq = 0.003 cm“^) and optically thick 
{fisd 3> 1), so that a great number interactions occur before the photon is lost. To simplify 
the system, we assumed that the refractive index does not change, we can visualize as if 
we were immersed in the environment. 

To encourage the development of flight we move the point of initial interaction ’’photons” 
to the center of the sample and we count N interactions from the origin to obtain the 
components a^, tty and of the resultant a (see Fig. 6). As discussed earlier, for a large 
number of interactions, applying the central limit theorem, these components must follow 
Gaussian statistics. 

A. Uniform System 

We study hrst the uniform system. In this case, the PDF for movement is given by equa¬ 
tion (6). Figure]^ shows histograms of the components resulting in, after 10 interactions. 
The vertical bars represent the histograms of displacement, and red curves, Gaussian 
functions that adjust data. We can see that although there are slight changes in the 
heights of the curves all have the same width oxo = 22.5 m. Glearly the components 
of the displacements resulting in good agreement with the expected Gaussian distribu¬ 
tion, which is consistent with the central limit theorem. This, despite the fact that we 
considered only 10 interactions. 


TABLE I. The system parameters which were performed the Monte Carlo calculations. 


layer 

n[-] 

Pa [cm 

fis [cm 

9[-] 

d [cm] 

sup. 

1.4 





optical system 

1.4 

0.003 

1867 

0.4 

10 

inf. 

1.4 
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FIG. 6. Histogram of the components ax, ay, and a^ of the resultant of the random walks for a 
uniform system. 
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Figure shows the histogram of the magnitude, a, of the resultant. As expected, the 
result hts very well to a Rayleigh PDF. 


TABLE II. The system parameters which were performed the Monte Carlo calculations. 


Optical properties 

Flights 

p [cm 

Ct [cm^] 

Ps{s) 

cte. 

cte. 

(pf) exp [ - {p.t)s\ 

cte. 

deltas 


cte. 

exp. neg. 

Ws 

(1 + pssy 

delta 

cte. 

(pt)exp(-(pt)s) 

delta 

deltas 

+ bp,2e~^^^ 

delta 

exp. neg. 

Ps 

(1 + PsS^ 

exp. neg. 

cte. 

(pt)exp(-(pt)s) 

exp. neg. 

deltas 

apie~^^^ + bp2e~^^^ 

exp. neg. 

exp. neg. 

Ps 

(1 + PsS^ 
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B. System with two types of particles 


Consider now the case of the medium with two types of particles. FDP for displacement 
is given by equation (16). Figure shows histograms of the resultant components after 
10 interactions. As in previous cases, vertical bars represent histograms and red curves 
Gaussian functions adjusted in height. This hgure was generated by taking a = 0.1 and 
a = 0.01, so that pi = O.l/i^ (long steps with a low probability) and /i 2 = 1.1/is (short 
step with high probability). As the difference between the values of the coefficients of 
scattering is great, it has the possibility of abrupt fluctuations. 

Unlike the previous case. Figure shows that it does not have a good £t to Gaussian 
curves. This means that after 10 steps, the statistics do not converge to such statistics. 
Should be noted that, if we increase the number of interactions or flights eventually 
expected convergence is obtained based on the central limit theorem. On the other hand, 
keeping the number of flights 10, but considering larger values of a (which implies that the 
two deltas distribution of the scattering coefficients are closer) is also obtained convergence 
Gaussian statistics. 


C. System with a negative exponential distribution of particles 


Now we assume that PDF governing displacement is given by equation (19). Figure 10 


shows the histogram of the components of the resultant, a, after 10 flights. As in previous 
figures, vertical bars represent histograms and red curves Gaussian functions adjusted in 
height .We see that the histograms of the components do not £t Gaussian curves.We can 
also see that the range of values ??taken by these components is much broader than in 
the previous cases, indicating that the fluctuations are much larger in flight, and can be 
up to an order of magnitude larger. 


This is not surprising, then flight statistics given by equation (19) represent a statistical 


approach to type Levy, and Levy type processes are caracterizazdos by violent fluctuations 
that make the resulting not converge to Gaussian statistics. 

We studied, however, the possibility of convergence after a very large number of interac- 
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FIG. 7. Histogram of the magnitude of the resultant a of the random walk. 
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FIG. 8. Histogram of the components ax, ay, and az of the resultant of the random walk, 


considering the PDF given by equation (15) with a 


0.1 and a = 0.01. 
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Ay, 


considering the PDF given by equation (19). 


and az of the resultant of the random walk, 


tions. Figure [TO] shows the histogram of the resultant components after 1, 000 flights. We 
see that, after such a large number of flights, the statistics converge if Gaussian statistics 
appear, although we should mention that this does not necessarily mean that the central 
limit theorem is valid in this type of situation. The statistics themselves seem to converge 
to Gaussian statistics, although we should mention that this does not necessarily mean 
that the central limit theorem is valid in this type of situation. 

The above results show that for this medium, if they occur a sufficiently large number 
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FIG. 10. Histogram of the components a^, ay, and the resultant of random walks for 1,000 
flights. 


of interactions, it will appear to be Gaussian statistics. However, such large fluctuations 
have important implications hnite system, as in hlms, in which the number of interactions 
is limited by the him thickness. 


V. COMMENTS AND CONCLUSIONS 

We have seen that the FDP adopted for the hights can determine the convergence or lack 
of it to Gaussian statistics. 

The fact of having non-Gaussian statistics and the ability to take big steps involve major 
changes in the properties of a him of this type of media. In these superdifusivos media, 
for example, the opacity of a him can be lowered considerably. 

The results show that in uniform media, after 10 steps, has a good convergence to Gaussian 
statistics. In the other two types of media considered more interactions are required to 
have these statistics. In particular, for the medium with negative exponential huctuations 
are required in the order of 1,000 interactions to approach these statistics. 
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